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Recently PAMELA collaboration published the cosmic nuclei and electron spectra 
with high precision, together with the cosmic antiproton data updated, and the 
Fermi-LAT collaboration also updated the measurement of the total e^e~ spectrum 
to lower energies. In this paper we develop a Markov Chain Monte Carlo (MCMC) 
package CosRayMC, based on the GALPROP cosmic ray propagation model to study 
the implications of these new data. It is found that if only the background electrons 
and secondary positrons are considered, the fit is very bad with x^cd ^ 3.68. Taking 
into account the extra e'^e~ sources of pulsars or dark matter annihilation we can 
give much better fit to these data, with the minimum x^cd ~ 0-83. This means the 
extra sources are necessary with a very high significance in order to fit the data. 
However, the data show little difference between pulsar and dark matter scenarios. 
Both the background and extra source parameters are well constrained with this 
MCMC method. Including the antiproton data, we further constrain the branching 
ratio of dark matter annihilation into quarks Bg < 0.5% at 2a confidence level. The 
possible systematical uncertainties of the present study are discussed. 

PACS numbers: 95.35.+d,96.50.S- 
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INTRODUCTION 



A very interesting progress made in the recent years in cosmic ray (CR) physics is the 
discovery of the excesses of positrons and electrons by several space- and ground-based ex- 
periments [l-5|. The positron and electron excesses challenge the traditional understanding 
of CR background. Many theoretical models were pro posed to exp 



ena, including the astrophysical scenarios (e 



dark matter (DM) scenario (e.g., 13l-ll6|) 



ain these new phenom- 
and see 12| for a review) and the 



Most recently Fermi-LAT team reported the updated measurements of the total spec- 
trum of electrons and positrons, extending to energies as low as several GeV 17| . PAMELA 
collaboration updated the observation of the p/p ratio and the absolute antiproton flux 18| . 
and reported the measurement of the pure electron spectrum at the first time 19|. The 
antiproton data extend to 180 GeV, without any hint of deviation from the background 
contribution [l^. For the PAMELA electron data, although there is no significant spectral 
feature above 30 GeV other than a single power-law, it is also consistent with models includ- 
ing extra e^e~ sources to explain the positron excess 19|. There is also hint of hardening of 
the electron spectrum compared with the low energy part (< 20 GeV), although the solar 
modulation may be important at these low energies. 

Since the accumulation of high-quality CR data, it is now important to extract more 
information from these data, i.e., estimating the CR background and the possible extra 
source parameters. Previously, one always constrained one or two parameters with other 
parameters fixed. This may lead to biased results, especially, when the parameters are 
strongly correlated. The global fitting procedure searches the maximum likelihood in the 
multiple dimensional parameter space rather than a reduced one and can give the posterior 
distribution by marginalizing other parameters in the Bayesian approach. The Markov Chain 
Monte Carlo (MCMC) procedure, whose computational time scales approximately linearly 
with the number of parameters, makes it possible to survey in a very large parameter space 
with the least computational cost. 

In our previous studies 20|, l2l|], we employed the MCMC method to fit the parameters of 
the DM scenario as well as the background parameters proposed to explain the e^e~ excesses. 
However, the propagation of CRs is treated with a semi-analytical way following Refs. 22 - 
24l |. A more precise description of the CR propagation is given by the numerical models, 



3 



such as GALPROP 



25| and DRAGON 26|, in which most of the relevant physical processes 



are taken into account, and the realistic astrophysical inputs like the interstellar medium 
(ISM) and the interstellar radiation field (ISRF) are adopted. There are many parameters in 
the CR propagation model, and it is very difficult to have a full and systematical survey of 
the parameter space. After embedding the numerical CR propagation tool into the MCMC 
sampler, we can use it to constrain the model parameters in a more efficient way. Recently 
';here are also several other works using the MCMC method to study the CR propagation 



27 



29|. 



In this work, we develo p a p ackage CosRayMC (Cosmic Ray MCMC), which is comprised 



of GALPROP, PYTHIA 



30| and MCMC sampler, to revisit the models to explain the 



positron and electron excesses and derive the constraints on the model parameters with the 
latest data. Two kinds of the extra e^e~ sources are considered: the pulsar scenario and 
the DM annihilation scenario. In both scenarios we consider the continuous distribution 
of the sources, although it is possible that one or several nearby pulsars or DM subhalos 
may explain the excesses 0, 0, lail-Q- As pointed out in Q, it was unlikely that the flux 



from any single pulsar was signiflcantly larger than that from others given a large number of 
known nearby and energetic pulsars. For DM subhalos the location and mass are uncertain, 
which also makes the constraints of the model parameters difficult. 

This paper is organized as follows. We briefly introduce the propagation of Galactic CRs 
in Sec. [Tll The description of the methodology is given in Sec. Illli The fltting results of 
the pulsar and DM scenarios are presented in Sec. IIVI Finally we give the conclusion and 
discussion in Sec. |Vl 



II. PROPAGATION OF GALACTIC COSMIC RAYS 

The charged particles propagate diffusively in the Galaxy due to the scattering with 
random magnetic fleld. There are interactions between the CR particles and the ISM and/or 
the ISRF, which will lead to fragmentation, catastrophic or continuous energy losses of the 
particles. For unstable nuclei the radioactive decay also needs to be taken into account. 
In addition, the overall convection driven by the stellar wind and reacceleration due to 
the interstellar shock will also affect the distribution function of CRs. For each species of 
particles we have a partial differential equation to describe the propagation process, with 
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the general form 



35| 



^ = g(x,p) + V-(D..V^-VeV^) + |-/D,,|-l^-|- U - |(V ■ Ve^)] (1) 

at op op p'' op I 6 J r/ 

where ip is the density of cosmic ray particles per unit momentum interval, Q{'x,p) is the 
source term, D^x is the spatial diffusion coefficient, Vc is the convection velocity, Dpp is 
the diffusion coefficient in momentum space used to describe the reacceleration process, p = 
dp/dt is the momentum loss rate, rj and r,. are time scales for fragmentation and radioactive 
decay respectively. Solving the partially coupled equations for all kinds of particles, we can 
get the propagated results of the CR spectra and spatial distributions. For more details 
about the terms in Eq. ([T]) please refer to the recent review paper 35 1. 

The secondary-to-primary ratios such as B/C and (Sc-|-Ti+V)/Fe, and the unstable- 
to-stable ratios of secondary particles such as ^°Be/^Be and ^^Al/^'^Al, are often used to 
constrain the propagation parameters because the ratios can effectively avoid the source 
parameters. Then one can use the spectra of the primary particles to derive the source 
parameters. There are some studies to constrain the propagation parameters based on the 



currently available data (e.g., [25| 



36 



37(). In 2^291 the MCMC method was adopted to 



fit both the propagation and source parameters of CRs. However, due to the quality of the 
observational data, the constraints on the propagation parameters are not very effective, 
and there may be also large systematic uncertainties. 

Since we focus on the electron and positron data in this work, which are not very effective 
to constrain the propagation parameters, we adopt the best fitting results of the propagation 
parameters given in 29[. The main propagation parameters are: -Do = 6.59 x 10^® cm^s~^, 
S = 0.30, va = 39.2 km s~^, = 3.9 kpc. The injection spectra of nuclei are adopted as 
7" = 1.91 (below 10 GV) and 72 = 2.40 (above 10 GV). What we should keep in mind is 
that there might be systematic errors of the determination of the propagation parameters 
(see e.g., jsSi), and the quantitative results of this work might be affected. 

We note that the newest measurements of the CR proton and Helium spectra by CREAM 



39| and PAMELA 40|] showed remarkable deviation of single power-law spectra. Further- 



more the spectral indices of proton and Helium are different. Such new features challenge the 



traditional understanding of the CR origin, acceleration and propagation (see e.g. j4l|, |42| 
for some explanations). In the present work we keep the study in the traditional frame and 
do not try to reproduce the detailed structures of the data. We find that the above adopted 
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injection parameters with proper adjusting of the normahzation can basically reproduce the 
PAMELA-CREAM proton data, as shown in the left panel of Fig. [H However, when com- 
paring with the Helium data, the model expectation seems to be systematically lower. We 
increase the relative abundance of Helium in GALPROP by 30%, and give roughly consis- 
tent result with PAMELA data (right panel of Fig. [1]). The CREAM data can still not be 
reproduced. Since the contribution to the secondary particles from CR Helium is only a 
small fraction, we do not think a higher Helium flux within a factor of 2 will substantially 
affect our following study. 




10' 10' 10 
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10' 10 10 
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FIG. 1: Proton (left) and Helium (right) spectra of the GALPROP model calculation, compared 
with the observational data. In each panel the higher line labelled "LIS" represents the local 
interstellar flux, and the lower one is the flux after the solar modulation. References of the data 
AMS 



are: proton 
iJ], CAPRICE94 



43[, BESS J4J], ATIC2 [45(L, PAMELA 
46|], ATIC2 y, PAMELA |4q], CREAM 



;4fl], CREAM 
39|. 



39j; Helium — BESS 



III. METHODOLOGY 
A. CosRayMC 

The CosRayMC (Cosmic Ray MCMC) code is built up by embedding the CR propagation 
code GALPROP (GALPROP v51) into the MCMC sampling scheme. The MCMC technique 
is widely applied to give multi-dimensional parameter constraints from observational data. 
Following Bayes' Theorem, the posterior probability of a model (which we refer to a series 
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of parameters 9) given the data (which are described by D) is 

v{e\D) ^mne), (2) 

where £i{9) = P{D\6) is the hkehhood function of the model 6 for the data D, and V{6) is the 
prior probabihty of the model parameters. In performing CosRayMC, MCMC is employed 
to generate a random sample from the posterior distribution V{6\D) which are fair samples 
of the likelihood surface. Based on the sample, we can get the estimate of the mean values, 
the variance as well as the confidence level of the model parameters. For details please refer 
to 0- 

The key improvement, compared with our previous study, is that the calculation of the 
likelihood C{6) are given by calling GALPROP to simulate the mock observations. By doing 
so, we can have a more precise description of the CR propagation. And as the better data- 
set is provided in the future, our code can also be used to make a determination of the CR 
propagation parameters. 

For the part with DM contribution to the CRs, the PYTHIA (PYTHIA v6.4) simula- 
;ion code is employed to calculate the final spectra of electrons, positrons and antiprotons 
3o| . Such spectra are then injected into the Galaxy and propagated with GALPROP. The 
PYTHIA code is also embedded in the CosRayMC code. 

B. Parameters 

We assume the injection spectrum of the background electrons to be a broken power-law 
function with spectral indices 71/72 below/above E'br- Note that for the shock acceleration 



scenario, the injection spectrum of particles can not be too hard 48|]. We set the priors 
that 7i > 1.5 and 72 > 1.5 in the MCMC scanning. The normalization of background 
electrons Abkg, taken as the flux of electrons at 25 GeV, is also regarded as a free parameter. 
For the background positrons and antiprotons, we adopt the GALPROP model predicted 



29| . Considering 



results with the best fltting source and propagation parameters given in 
the fact that there are uncertainties about the ISM density distribution, the hadronic inter- 
action model, and the propagation parameters determined from the secondary-to-primary 
ratio data, we will further employ factors Ce+ and Cp to rescale the absolute fluxes of these 
secondary particle. 
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For energies below ~ 30 GeV the solar modulation effect is important and needs to 
be considered. In this work we adopt the force-field approximation to calculate the solar 
modulation 49|. The modulation potential depends on the solar activity. For the period 



which PAMELA works the modulation potential is estimated to be = 450 — 550 MV |40 |. 
In our MCMC fit, the modulation potential is also taken as a free parameter. Thus, for 
the background model we have 7 parameters in total 



"^bkg = {7l, 72, -Ebr, Akg, 4>, Ce+, Cp}. 



(3) 



Pulsars are thought to be the most natural candidates to generate high energy positrons 



and electrons through the cascade of electrons accelerated in the magnetosphere [31|, [50 1. 
The spectrum of e~^e~ escaped from the pulsars can be parameterized as power-law with a 
cutoff at Ec, dN/dE oc ApsrE~°' exp{—E/ E^), where the power-law index a ranges from 1 



to 2.2 according to the radio and gamma-ray observations 



311 ]. The cutoff energy of the 



injected e~^e ranges from several 



and parameters of the pulsars 
following form 51] 



34 



tens GeV to higher than TeV, depending on the models 



f{R,z)cx 



B_ 

'Rc. 



50| . For the spatial distribution of pulsars, we adopt the 



■0 



exp 



Rc. 



'0 



exp — 



(4) 



where Rq = 8.5 kpc is the distance of solar system from the Galactic center, Zg ~ 0.2 kpc 
is the scale height of the pulsar distribution, a = 2.35 and b = 5.56. There will be a further 
normalization factor Apsi- 

Alternatively, DM annihilation or decay models are widely employed to explain the e'^e~ 
excesses. Since to date one can not distinguish DM annihilation from decay with the e~^e~ 
data 



20 



211 ]. we only consider the annihilation scenario here. Considering the non-excess 



of 



15 



AMELA antiproton data, the DM annihilation final states need to be lepton-dominated 



16| . Similar as the previous works 20|, l2l|, we do not employ the annihilation final 



states based on any DM models, but to constrain the final states from the data we use a 
model-independent way instead. Such constraints will be helpful for the understanding of 
the DM particle nature. 

The annihilation final states are assumed to be two-body e^e^ , fi^fi^ , t^t~ and gg^, with 



^ Since the antiproton productions of various quark flavors do not differ much from each other, here we do 
not distinguish quark flavors but adopt uu channel as a typical one of quark final states. 
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branching ratios B^,, B^, Bj- and Bq respectively. It was also shown that the interactions 
with low-mass intermediate bosons which then decay to lepton pairs could fit the CR e^e" 
data well 



52 



53| . However, since the spectral shapes of the electrons and positrons from 



these models do not have distinct properties which can be easily distinguished from fi~^fi~ 
and r+r^, we do not include such final states in this study. The resulting positron, electron 
and antiproton spectra of these final states are calculated using the PYTHIA simulation 



package 



30|. 



The DM density profile is assumed to be an Einasto type 54 1 

p(r)=P-.exp[-^(^-l)], (5) 

where a ~ 0.17, r„2 ~ 15.6 kpc and p_2 ~ 0.14 GeV cm~'^ according to the recent high 
resolution simulation Aquarius |55|. The corresponding local density of DM is about 0.44 
GeV cm~^, which is consistent with the results of a larger local density derived in recent 



studies 



56N58|. Then the source function of positron produced from DM annihilation is 

av) 



q{E,r) 



2m?. 
^ f 



X P W, 



(6) 



where is the mass of the DM particle, {av) is the velocity weighted annihilation cross 
section, ^ | ^ is the positron production rate from channel / of one annihilation. 
The most general parameter space is summerized as below. 



tol 



{"Pbkg}, background 
{'Pbkg,^psr,a,^c}, pulsar 
{■Pbkg, Cj5, m^, {av),Be, fi^, Br, Bu}, DM. 



(7) 



For the DM annihilation scenario, we further consider two cases. One is to use the e~^e~ 
data only in the fit, which can be directly compared with the pulsar scenario. In the other 
case, the antiproton data are included in order to constrain the coupling to quarks of DM. 
In the case without the antiproton data, Bu and Cp are set to zero. 
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IV. RESULTS 



A. Background 



First we run a fit for tlie pure background contribution. The data included in tlie fit 
are PAMELA positron fraction PAMELA electron 19|, Fermi-LAT total e+e" 17|, and 



HESS total e+e- 



\i 



For the PAMELA positron fraction data, we only employ the data points with E > 5 GeV 
in the fit. For the lower energy part we can see that the PAMELA data are inconsistent with 
the previous measurements. This may be due to more complicated solar modulation effect 



or even the charge dependent modulation 



59 



611 ]. It was shown that either a model with 



different modulation potentials for positive and negative charged particle 



601 ] or a detailed 



Monte Carlo approach to solve the stochastic differential equations of the e~^e motion 



6l| can give consistent description to the PAMELA and previous data. The demodulated 
interstellar spectra of electrons and positrons are also consistent with the conventional CR 
background model expectation. Here we do not consider the detailed solar modulation 
models. But note that the solar modulation model may affect the quantitative fitting results. 
This is one kind of systematical errors. 

The best-fit results of the positron fraction and electron (or e~^e~) spectra are shown 
in Fig. [21 The best-fit parameters are compiled in Table [H with the minimum Xmd ~ 
415.4/113 = 3.68. Such a large reduced means that the fit is far from acceptable but 
systematics dominated. Using a goodness-of-fit test we find that the background only 
case is rejected with a very high significance ~ 20a. That is to say the data strongly favor 
the existence of additional degrees of freedom. From Fig. [2] we can clearly see that the 
background model under-estimates both the positron fraction and the total e~^e~ fluxes. 
The extra sources of e~^e~ are necessary to explain simultaneously the positron fraction and 
total e~^e~ data. 

The poor fit makes nonsense to discuss the physical implication of the parameters. There- 
fore we leave the discussion of the parameters in the next subsections, when the extra sources 
of e^e~ are taken into account. 
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FIG. 2: Best-fit positron fraction (left) and electron spectra (right) of the pure background case. 
Note that in the right panel the PAMELA data are for pure electrons, while other data are for 
the sum of electrons and positrons. References of the observational data are: positron fraction — 




£21, HEAT94+95 
Fermi-LAT 



63], HEATOO 



Mi, PAMELA electron ~ PAMELA pj], ATIC 



B. Pulsar scenario 



Compared with the conventional GALPROP background model, 71 = 1.63 ± 0.08 is 
consistent with 1.60 as given in 65|. Note that there is a correlation between parameters 

71 and the modulation potential 0, as shown in Fig. |3l The high energy injection spectrum 

72 = 2.74 ± 0.02 is softer than 2.54 as adopted in the conventional model 65|]. In the 
conventional model only the background contribution is employed to fit the data, while 
here an additional component from the extra sources are added together. Thus the high 
energy spectrum can be much softer. This might be important for the understanding of the 
background contribution to the CR electrons. 

There is a factor Ce+ ~ 1.6 needed for the background positrons to fit the data. As 
discussed above, such a factor may be ascribed to the uncertainties of the propagation 
parameters, the ISM density and the strong interaction cross sections. Those uncertainties 
may be energy-dependent and can not be simply described with a constant factor (see e.g., 
jsj). For example in this work we use the parameterization given in 66 1 (Kamae06) to 
calculate th e interaction to generate positrons. Compared with the pp collision model 
Badhwar77 67| as adopted in GALPROP, the Kamae06 model gives systematically fewer 
positrons, especially for energies from several to tens of GeV jsj. If these input uncertainties 
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TABLE I: Fitting parameters with la uncertainties or 2a limits. Note that for the "bkg" case 
the reduced is too large that the uncertainties of the parameters should not be statistically 
meaningful. 





bkg 


bkg±pulsar 


bkg±DM (without p) bkg±DM (with p) 


71 


< 1.535(95%C.L.) 


l.UOU_Q Qg3 


1 614+°-°^2 

—0.070 


< 1.610(95%C.L.) 


72 


2.559 ±0.008 


2.738 ±0.021 


2.730 ±0.018 


2.706 ±0.013 


log(Akg) ^ 


-8.958 ± 0.004 


-8.981 ±0.010 


-8.984 ± 0.008 


-8.997 ± 0.006 




o t;c;4+0.133 
<J.OO4_0.i29 


j.yDy_o.29i 


3.99li°:i° 


4.283l°J^6 


cPiGV) 


0.333 ±0.019 


0.499 ± 0.071 




0.371 ±0.037 


Ce+ 


1.473 ±0.037 


1.596 ±0.111 


1-516^°:°8^^ 


1.394 ±0.053 


Cp 


— 


— 




1.210 ±0.045 


log(Apsr) ^ 


— 


-28.36l°j3 






a 


— 








Ec 




n 707+0.154 
'^•'^'-0.164 






M^{TeV) 






2.370l°:i^ 


9 o/i-i +0.492 


log[av{cm^ s~^)] 






-22.33 ±0.13 


-22.34 ±0.13 


Be 






< 0.369(95%C.L.) 


< 0.379(95%aL.) 








< 0.381(95%C7.L.) 


< 0.334(95%C7.L.) 


Br 






0.688i°:lf, 


Q 71 O+0.141 

' ■'-"J-0.152 


Bu 








< 0.005(95%C7.L.) 


xVd.o.f 


3.676 


0.833 


0.827 


1.078 



^ ^bkg is in unit of cm "^sr ^MeV 
^ A-psr is in unit of cra^"^ sr^^ M eV~ 



are better understood in the future, we can include the background positron production 
with a more general form in the MCMC fit. 

The fitted parameters for pulsars are a = 1.201]^q;}295 = 0.797^0.164 TeV. The con- 
straints of the pulsar parameters are weaker when compared with the background param- 
eters. This is because the cutoff energy is mainly determined by HESS data. However, 
for HESS data, the very large systematic errors due to poor absolute energy calibration 
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FIG. 3: Two dimensional constraints of part of the model parameters for the pulsar scenario. 
In off-diagonal figures, the inner contour denotes the 68% confidence level (C.L.) while the outer 
contour stands for 95% C.L. The diagonal figures are one dimensional probability distributions for 
corresponding parameters. 



(not shown in Fig. H]) makes it impossible to precisely determine E^. In the calculation 
the relative systematic uncertainty of the flux is estimated to be (F — 1)AE/E [l^]. For 
~ 15% uncertainty of the energy scale, the flux uncertainty is ~ 30% and ~ 45% for energies 
below and above ~ 1 TeV. Due to the correlations between E^ and a, Ap^^, the other two 
parameters are also relatively less constrained. 
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FIG. 4: Same as Fig. [Hbut for the model with pulsars as the extra sources of e+e . 



C. Dark matter annihilation scenario 

First we consider the case without including the antiproton data. We found that the 
background parameters are similar with the pulsar scenario. The overall value of the 
best-fit is about 90.2, which is very close to Xpsr ~ 91.6. That is to say both the astrophysical 
scenario and DM scenario can give comparable fit to the e^e~ data. We may not be able 
to discriminate these two scenarios from the e^e~ spectra 



other probes such as the e^e anisotropy [72| and 7-rays [73|. 



34 



68 



71| . but need to resort to 



The fitted mass and cross section of DM are = 2.37lo;39 TeV, (av) = 4.7+J:2 x lO'^^ 
cm^ s"^. Similar as the pulsar case, these parameters are also less constrained due to 
the large uncertainties of the HESS data. We can also determine the branching ratios to 
different flavors of leptons: B^{2a) < 0.37, B^X'^a) < 0.38 and = OmtoH These results 



20|, but with larger uncertainties because 



are qualitatively consistent with our previous fit in 
we include the systematic errors of HESS data in the present work. We should be cautious 
that the determination of the branching ratios may suffer from systematic uncertainties 
of the background electrons and positrons. As discussed in the previous sub-section, the 
present understanding of the background positrons is far from precise. Therefore the quoted 
fitting uncertainties of the branching ratios may be under-estimated. 

Then we consider the case with antiproton data. The coupling to quarks of DM an- 
nihilation is taken into account. The two dimensional constraints of part of parameters 
are present in Fig. [5l Compared with the case without antiproton data, we find that the 
background parameters are slightly different. The reason for the change is mainly due to 
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FIG. 5: Same as Fig. [3] but for the DM scenario with antiproton data. 



the low energy spectrum of the background antiprotons. As shown in Fig. [6] the calculated 
secondary (including tertiary) antiprotons are basically consistent with the PAMELA data 
for energies higher than several GeV, but under-estimate the flux for the low energy part. 



This problem was pointed out years ago 



78|, that in the diffusive reacceleration model 



which gives the best fit to the B/C data the antiprotons are under-estimated. There seems 
to be a contradiction between the B/C and the antiproton data^. To better understand 



^ Note, however, in 
to the B/C data 



791 

36| 



the calculated antiproton flux based on the propagation parameters fitted according 
was consistent with the observational data. As pointed out in [78[ the fit in [36| 
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FIG. 6: Fitted 2a range of the antiproton 
the data are: AMS 



7|, BESSOO 
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fluxes of the DM annihilation scenario. References of 
BESS02 



Yg], PAMELA 



1^. 



this issue we may need more precise measurement about the B/C data. Back to this work, 
a lower antiproton flux at low energies will require a smaller solar modulation potential, 
and therefore the background parameters 71, 72 and E'br will change accordingly due to the 
correlations among them as shown in Fig. [5l 

One thing important is that we can derive a self-consistent constraint on the quark 
branching ratio of DM annihilation. The marginalized 2(7 upper limit of Bu is about 0.5%. 



Compared with our previous study 



21| the upper limit of i?„ is several times smaller. This 



is probably due to the bad fit to the antiproton data with the background model. We have 
run a test to employ the background antiproton spectrum used in 2l|] and found that the 2a 
upper limit of Bu is about 2.1%, which is consistent with the previous results. This means 
that the current constraint on Bu is systematics dominated. Better understanding of the 
background contribution to the antiproton flux is necessary to further address this issue. 
We plot the 2cr range of the antiproton fluxes in Fig. [61 



was actually based on the high energy data. And in their results the low energy B/C was in fact over 
estimated. Furthermore the antiproton production cross sections would also lead to uncertainties [78 1. 



Alternatively, in [80[ a unified model to explain the B/C and antiproton data was proposed, with an 



empirical modification of the diffusion coefficient at low speed of particles. 
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V. CONCLUSION AND DISCUSSION 

Recently more and more observational data of CRs with unprecedented precision are 
available, which makes it possible to better approach the understanding of the basic problems 
of CRs. Based on MCMC analysis in our previous work 
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21|, we embed GALPROP and 



PYTHIA into MCMC sampler and study the implication of the newest CR data, including 
the positron fraction, electrons (pure e~ and e~^e~) and antiprotons from PAMELA, Fermi- 
LAT and HESS experiments in this work. 

We work in the frame of diffusive reacceleration propagation model of CRs. The prop- 



29|, 
40|. 



agation parameters are adopted according to the fit to currently available B/C data 
with a slight adjustment of the Helium abundance to better match the PAMELA data 
We find that the pure background to explain the CR e+e~ data is disfavored with a very 
high significance. Therefore it is strongly implied that we may need some extra sources to 
produce the positrons/electrons. 

We then consider two different scenarios, which are widely discussed in recent literature, 
to explain the e+e~ excesses. One is the astrophysical scenario with pulsars as the typical 
example, and the other is the DM annihilation scenario. Using the global fitting method, 
we can determine the parameters of both the background and the extra sources. The low 
energy spectral index of the background electrons is consistent with the conventional model 
adopted in the previous study, while the high energy index is softer in this work. We find 
that, together with the background contribution, both of these scenarios can give very good 
fit to the data. For the case with only the e~^e~ data, the goodness-of-fit of the pulsar 
scenario and DM scenario are almost identical. With the antiproton data included, we can 
further constrain the coupling to quarks of the DM model. The branching ratio to quarks 
of the DM annihilation final states is constrained to be < 0.5% at 2cr confidence level. 

The constraint on the quark branching ratio (or the cross section to quarks) allows us 
to make a comparison with the direct detection experiments, assuming some effective in- 



teraction operators between DM and standard model (SM) particles. As shown in 81| . 
generally the constraint of DM-SM coupling from direct experiments is much stronger than 
the indirect search for spin-independent interactions. However, for the spin-dependent in- 
teractions, the indirect search constraint could be comparable or stronger than the direct 
experiments. In Fig. [7] we give the constraints on the spin-dependent DM- neutron scattering 
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cross section according to the cross section of DM annihilation to quarks, assuming the axial 



vector interaction form of DM particles and SM fermions. The ratio between {av) and a. 



SD 
xn 

is 5c(l + ran/nix)'^ /m^^^ y^l — m'^/m'^mg [8l|, where c is light speed, m„ is the mass of 
neutron and rriq is the mass of quark. The sum is for all the flavors of quarks. Considering 
the uncertainties of the background antiprotons, we show the results of two antiproton back- 
ground models with shaded region: the result calculated in this work and the one used in 



2l| . The results show that the constraint on the spin-dependent cross section between DM 



and nucleon from PAMELA antiproton data is much stronger than the direct experiments. 



10 



10 



-36 



-37 



E 



10 



-40 



10" 



10 



PAMELA pbar 
CDMS-2006 
XEN0N1 0-2008 



10' 10^ 
m (GeV) 



10 



FIG. 7: Constraints on the spin-dependent DM- neutron scattering cross section from the PAMELA 
antiproton data assuming axial vector interaction form of DM particles and SM fermions. The 
shaded region represents the uncertainties of the astrophvsical background of antiprotons. Also 



shown are the results from CDMS 



8^ and XENONIO 



Although the existence of the extra sources is evident, we still can not discriminate the 
pulsar and DM scenarios right now. These two models both can fit the data very well, with 
Xred — 0.833 and Xrcd — 0.827 respectively. Other new probes such as the e~^e~ anisotropy 
and 7-rays are needed to distinguish these models. 

The 7-rays should give further constraints on the models, especially for the DM annihi- 
lation scenario. However, the diffuse 7-ray data of Fermi-LAT are still in processing, and 
the modeling of the 7-ray data seems not trivial. Furthermore the 7-rays are more sensitive 
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to the density profile of DM, wliich is not tlie most relevant parameter of the present study. 
Therefore we do not include the 7-ray data. It should be a future direction to include the 
7-ray data in the CosRayMC package. 

There are some systematical uncertainties of the current study, which are mainly due to 
the lack of knowledge about the related issues and can be improved in the future. First, 
the solar modulation which affects the low energy spectra of CR particl es [ E < 30 GeV) 
if not clear. In this work we use the simple force-field approximation 49| to deal with 
the solar modulation. But such a model seems to fail to explain the low energy data of 
;he PAMELA positron fraction, which may indicate the charge dependent modulation effect 



59 



61| . Second, the propagation model is adopted as the diffusive reacceleration model with 



the parameters best fitting the B/C data 29|. There are uncertainties of the propagation 
parameters. Furthermore the diffusive reacceleration model may also have some systematical 
errors when comparing with all of the CR data. For example the antiprotons below several 
GeV is under-estimated in this model. We need to understand the propagation model better 
with more precise data. Third, the secondary positron and antiproton production may suffer 
from the uncertainties from the ISM distribution and the hadronic cross sections. All of 
these uncertainties may affect the quantitative results of the current study. Nevertheless, 
with more and more precise data available in the near future (e.g., AMS02, which was 
launched recently), and better control of the above mentioned systematical errors, it is a 
right direction to do the global fit to derive the constraints and implication on the models. 
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